IDW Subroutine

SUBROUTINE IDW &
!
(network, grid, method, power, near)


IMPLICIT NONE

!Arguments with intent (in):
TYPE (ObservationalNetwork), INTENT(IN) :: network
INTEGER (KIND = short), INTENT (IN) :: method
REAL (KIND = float), OPTIONAL, INTENT (IN) :: power
INTEGER (KIND = short), OPTIONAL, INTENT (IN) :: near

!Arguments with intent (inout):
TYPE (grid_real), INTENT (INOUT) :: grid

!Local declarations
TYPE (ObservationalNetwork) :: activeNetwork
INTEGER (KIND = short) :: count
INTEGER (KIND = short) :: i, j, k, s, t
INTEGER (KIND = short) :: nearPoints
REAL (KIND = float)    :: actPower
REAL (KIND = float),POINTER   :: dist(:,:)  
REAL (KIND = float)   :: distIJ
REAL (KIND = float)   :: denom 
REAL (KIND = float)   :: weight 
!-----------------------end of declarations------------------------------------

!check for available measurements
CALL ActualObservations (network, count, activeNetwork)

!allocate distances vector
ALLOCATE (dist (activeNetwork % countObs + 1,2) )

!set points
point1 % system = grid % grid_mapping
point2 % system = grid % grid_mapping

!set near
IF (PRESENT (near)) THEN
  IF (activeNetwork % countObs <= near) THEN 
    nearPoints = activeNetwork % countObs
  ELSE
    nearPoints = near
  END IF
ELSE
  nearPoints = activeNetwork % countObs !use all data
END IF

!set power
IF (PRESENT (power)) THEN
  actPower = power
ELSE
  actPower = 2. !default to square
END IF

DO i = 1, grid % idim
  DO j = 1, grid % jdim
    IF (grid % mat (i,j) /= grid % nodata) THEN
      dist = -9999.99
      !compute distance from cell center to observations
      CALL GetXY (i,j,grid, point1 % easting, point1 % northing)
      DO k = 1, activeNetwork % countObs        
         point2 % northing = activeNetwork % obs (k) % xyz % northing  
         point2 % easting = activeNetwork % obs (k) % xyz % easting
         distIJ = Distance (point1,point2)
                  
		 !search for neighbours
		 DO s = 1, nearPoints
			IF ( dist(s,1) == -9999.99 .OR. &
			     dist(s,1) > distIJ       ) THEN
				DO t = nearPoints, s, -1
					dist(t+1,1) = dist(t,1)
					dist(t+1,2) = dist(t,2)
				END DO
				dist(s,1) = distIJ
				dist(s,2) = activeNetwork % obs(k) % value
				EXIT
			END IF
		 END DO
	  END DO	 
	  !compute denominator
	  denom = 0.0
	  IF (method == 1) THEN !Shepard's method
	      DO s = 1, nearPoints
	 	      IF (dist(s,1) <= 1.) THEN
			     !observation in cell, set minimum distance to avoid dividing by zero
			     dist(s,1) = 1.
		      END IF
		      denom = denom + dist(s,1) ** (-actPower)
	      END DO
	      ! weighted value
	      grid % mat(i,j) = 0
	      DO s = 1, nearPoints
	        weight = dist(s,1) ** (-actPower) / denom
	        grid % mat(i,j) = grid % mat(i,j) + weight * dist(s,2)
	      END DO
	   ELSE !Franke & Nielson method
	      DO s = 1, nearPoints
	 	      IF (dist(s,1) <= 1.) THEN
			     !observation in cell, set minimum distance to avoid dividing by zero
			     dist(s,1) = 1.
		      END IF
		      denom = denom + ( (dist(nearPoints,1) - dist(s,1)) / (dist(nearPoints,1) * dist(s,1)) ) ** actPower
	      END DO
	      ! weighted value
	      grid % mat(i,j) = 0
	      DO s = 1, nearPoints
	        weight = ( (dist(nearPoints,1) - dist(s,1)) / (dist(nearPoints,1) * dist(s,1)) ) ** actPower / denom
	        grid % mat(i,j) = grid % mat(i,j) + weight * dist(s,2)
	      END DO
	   END IF		 
    END IF
  END DO
END DO

DEALLOCATE(dist)
DEALLOCATE(activeNetwork % obs)

RETURN

END SUBROUTINE IDW